Dark soliton dynamics in Bose-Einstein condensates at finite temperature 
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The dynamics of a dark soliton in an elongated Bose-Einstein condensate at finite temperatures 
is studied using numerical simulations. We find that in the presence of harmonic confinement the 
soliton may oscillate even at finite temperatures, but with an amplitude that increases with time, 
indicating the decay of the soliton. The timescale of this decay decreases both with increasing 
temperature and with increasing initial soliton velocity. Simulations performed for the experiment 
of S. Burger et at, Phys. Rev. Lett. 83, 5198 (1999), reveal excellent agreement with the observed 
soliton decay, confirming the crucial role of the thermal cloud in soliton dynamics. 
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Bose-Einstein condensates (BECs) in ultracold gases 
provide an ideal test ground for nonlinear physics. At 
temperatures close to absolute zero, the condensate 
can be accurately described by the Gross-Pitaevskii 
(GP) equation which has the form of a nonlinear 
Schrodinger equation where the nonlinear term arises 
from interactions between atoms. A similar equation ap- 
pears in nonlinear optics, where it is known to support 
soliton solutions 0. It is natural, therefore, to expect 
solitons in BECs. 

Both bright and dark solitons have indeed been ob- 
served in experiments. Bright solitons can be formed 
when the atomic interactions are attractive [H 0, 
whereas "gap" bright solitons have been created experi- 
mentally in repulsive condensates in optical lattices Q. 
Repulsive condensates in harmonic traps, however, can 
only support dark solitons, which correspond to prop- 
agating one-dimensional localized minima in the den- 
sity. Dark solitons have been generated experimen- 
tally 0, H, EJ EH j and their zero temperature properties 
have been studied theoretically in both three-dimensional 
geometries, where they exhibit dynamical instabilities 
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13, 14, 15, 1 



12j . as well as in one dimension 
Although modelling the condensates at zero temper- 
ature provides important insight into dark solitons, in 
practice experiments are always performed at non-zero 
temperatures, where the presence of a cloud of non- 
condensed particles provides a mechanism for the decay 
of the soliton. Despite the experimental evidence for such 
dissipative effects, the finite temperature properties of 
solitons have received relatively little attention. Refs. 
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18| predicted dissipation at finite temperatures by 
considering the reflection of excitations from a soliton in 
a uniform condensate. The role of quantum fluctuations 
has also been studied fl9j |. 

This Letter uses numerical simulations to perform a de- 
tailed quantitative study of the dissipative dynamics of 
dark solitons in elongated harmonic traps at finite tem- 
peratures. At zero temperature, a dark soliton is pre- 
dicted to oscillate in the axial direction [13, 13, 0, 0, 
[l6l |. However, the presence of a thermal cloud leads to 



damping. As a result, the depth of the propagating soli- 
ton decreases in time, leading to an increase in the am- 
plitude of the oscillations, which eventually approaches 
the half-length of the condensate. This is consistent with 
a decrease in the soliton energy as a function of time. We 
perform detailed investigations of this decay for different 
temperatures and initial soliton depths, and separately 
assess the role of mean field coupling and binary colli- 
sions between the atoms. In particular we directly sim- 
ulate the experiment of Burger et al. [7J , where the evo- 
lution of dark solitons, created by phase imprinting, was 
observed using time-of-flight absorption imaging. Our 
simulations show that the soliton effectively disappears 
after only half an oscillation (Fig. [IJ bottom images), 
in good agreement with the experimental findings, thus 
proving the crucial role of the thermal cloud. This should 
be contrasted with the undamped dynamics predicted by 
the GP equation (top images). 
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FIG. 1: Evolution of a dark soliton at T = (top) and 
T — 0.5 T c (bottom), corresponding to the experimental pa- 
rameters of Ref. [7J. In each simulation we initially imprint 
a phase on the condensate, then allow it to evolve for the 
time shown before releasing it from the trap. Plotted are 
condensate column densities after a subsequent expansion of 
t = 4 ms. 

Our simulations are based on the formalism of 
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Zaremba, Nikuni, and Griffin (ZNG) [20j, where the dy- 
namics of the condensate and the thermal cloud are de- 
scribed by the coupled equations: 
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Eq. |T]) is a generalized GP equation for the condensate 
wavefunction *£(r,t), while Eq. © is a Boltzmann equa- 
tion for the thermal cloud phase space density f(p,r,t) 
(where p is the momentum). The condensate and ther- 
mal cloud densities are defined as n c — \^\ 2 and n — 
J fd 3 p/h 3 respectively, while g — 4irh 2 a/m parameter- 
izes the mean field interactions between atoms of mass 
to and scattering length a. The effective potential is 
U = V + 2g(n c + n), with V — m(u\r 2 + w 2 z 2 )/2 repre- 
senting the harmonic trap. Apart from the mean field 
coupling between the condensate and non-condensate, 
Eq. fl} contains a term R — (h/2n c ) J C\i d 3 p/h 3 aris- 
ing from collisional processes, C%2, that add or remove 
a single particle from the condensate. The C22 term in 
Eq. |(5J) denotes binary collisions involving only thermal 
atoms. 

This theory has already been used to model the 
damping of collective modes, demonstrating good agree- 
ment with experiments (see e.g. Ref. plj and references 
therein) , with the numerical methods used to solve these 
equations described in Ref. [23 ]. Here, we also model the 
dynamics of the thermal cloud using N body simulations. 
Eq. §Q is solved in cylindrical coordinates (r, z) using 
a Crank-Nicholson scheme, with the singularity in the 
Laplacian at r = avoided by offsetting the radial grid 
by half of a grid spacing. 

For numerical convenience, our initial analysis is per- 
formed for a relatively small sample of N = 2 x 10 4 87 Rb 
atoms in a trap with frequencies lo z = 2~k x 10 Hz and 
uj±_ = 2ir x 2500 Hz. The cloud is thus highly elongated 
along the axial direction, with the tight transverse con- 
finement suppressing "snake-instabilities" which lead to 
decay of the soliton into vortices (UEUEl 12 1. 

The initial condition for each simulation is established 
by first self-consistcntly finding the equilibrium state of 
the condensate wavefunction and thermal cloud Bose dis- 
tribution f(p, r) = {exp[(p 2 /2m + U - n)/k B T] - l}" 1 
(where fi is the condensate chemical potential) at a given 
temperature. A soliton is then generated by multiply- 
ing the ground state condensate wavefunction ^(r) with 
& a (z) = /3tanh(/3z/£) + i(v/c), where (3 = yjl - (v/c) 2 , 
£ = h/ \/mgn~ c is the condensate healing length, and 
c = y/ gn c j (2m) is the speed of sound. 

Firstly, we consider only the effect of mean-field cou- 
pling (n c , h) between the condensate and thermal cloud, 
setting R = C12 = C22 = in Eq. We perform sim- 
ulations for a range of initial velocities v, where smaller 



v is associated with deeper solitons (with a completely 
dark soliton being stationary). To study the soliton de- 
cay process in detail, we first consider a rather idealized 
example of a deep (slow) soliton with v = 0.1c. The sub- 
sequent dynamics is illustrated in Fig. [21 which shows the 
evolution of the longitudinal density profile for three dif- 
ferent scaled temperatures, T/T c , where T c is the critical 
temperature for BEC. The density inhomogeneity due to 
the harmonic trap is evident, with the soliton appearing 
as a dark line due to its low density. 

Fig. [2^a) illustrates the dynamics at T = 0, where 
the soliton oscillation has a constant amplitude and a 
frequency close to w z /\2, in agreement with the predic- 
tions of the GP equation [l2|j HI 14 1. In this case, the 
soliton maintains an almost constant depth, and the ex- 
trema of the oscillation occur when the density at the 
soliton minimum reaches zero. Fig- EKb) shows the result 
at T = 100 nK (0.2 T c ), where a small thermal cloud is 
also present. The steady increase in the oscillation am- 
plitude accompanies a decrease in the soliton depth. At 
T = 200 nK (0.4 T c ) [Fig. [2[c)] the rate of increase in 
amplitude is larger, and the soliton quickly becomes very 
shallow, making it difficult to visualize in the image. 




FIG. 2: (Color online) Temporal evolution of the condensate 
density as a function of z along a cross-section at r = Ar/2, 
where Ar is the grid spacing. Light colors represent high 
densities and dark low densities, with the soliton appear- 
ing as a dark line. Both the position and time units are 
expressed in terms of the axial trap frequency ui z , where 
a z = [h/ (mw z )f 2 . The three different plots are for (a) T = 0, 
(b) T ~ 0.2 T c , and (c) T ~ 0.4 T c , where T c ~ 500 nK, and 
the initial velocity is v — 0.1c. 



This increase in oscillation amplitude is a consequence 
of damping of the soliton due to dissipation. This is most 
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clearly seen by considering the energy of the soliton 
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evaluated at r — Ar/2, which is the nearest grid point 
to the z-axis. The integral is taken over a small interval 
28 centered on the soliton position, z s , where we choose 
5 = 5£. The background density for a condensate without 
a soliton is represented by n . 

In Fig. [3] we plot the soliton energy as a function of 
time. We see that for T = (dotted black line, top) the 
energy is constant, apart from an oscillation arising from 
continuous soliton-sound interactions (lBj : an accelerat- 
ing soliton in a harmonic trap periodically emits sound, 
leading to an instantaneous loss in energy. However, since 
the condensate is of finite size, the sound cannot escape 
from the system and is therefore periodically re-absorbed 
by the soliton, thus leading to no net decay in energy. In 
contrast, at finite temperatures there is a gradual loss of 
energy. While this effect is only marginal at relatively 
low temperatures (T ~ 0.2 T c ) represented by the black 
solid line, the dissipation is significantly enhanced with 
increasing temperature, as can be seen from the black 
dashed line (T ~ 0.4 T c ). 
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FIG. 3: Soliton energy, E s , as a function of time for the 
parameters of Fig. O and temperatures T — (dotted), T ~ 
0.2 T c (solid) and T ~ 0.4 T c (dashed). In the latter two cases, 
black lines plot the collisionless simulations, while gray lines 
include collisions. 

Although the collisionless regime discussed thus far 
captures some of the essential physics, experiments with 
solitons are actually performed in the presence of colli- 
sions (C\2 7^ 0, C22 7^ 0), to which we now turn our 
attention. The energy decay when collisions are included 
is represented by the gray lines in Fig. [3] Comparing to 
the collisionless results, one sees that the decay rate is 
greatly enhanced, and for T ~ 0.4 T c the energy falls to 
zero already at u z t ~ 20. 



Simulations for the soliton dynamics have been per- 
formed for a range of temperatures, initial velocities, and 
with and without collisions. In order to quantify the re- 
sulting soliton decay, we define a "half- life", 7i/ 2 , which 
corresponds to the time required for the soliton depth to 
reach half of its initial value. This quantity is relevant 
since it is closely related to the contrast of the absorp- 
tion images obtained experimentally after time-of-fiight 
expansion. The dependence of the half-life on the above 
parameters is shown in Fig. [3] The procedure for ob- 
taining T1/2 is illustrated in the inset, which plots the 
time-dependence of d, the depth of the soliton normal- 
ized to its initial value. Dashed lines show the point at 
which d = 0.5 for the first time, while the dotted lines in- 
dicate the points at which d = 0.45 and d — 0.55. These 
enable us to add error bars to our half-life estimates, in 
a manner which could account for experimental uncer- 
tainties in the contrast of the absorption images. Note 
that these error bars may be asymmetrical due to the 
non-monotonic change of the soliton depth arising from 
soliton-sound interactions. 
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FIG. 4: Time, Tim, for the soliton depth to decay to half of its 
initial value as a function of scaled temperature, T/T c , for dif- 
ferent initial velocities v = 0.1c (circles), v — 0.25c (triangles), 
and v — 0.5c (diamonds). Solid (open) symbols correspond to 
collisionless (collisional) simulations. The dashed line marks 
the period of one soliton oscillation, Inset: 
Time- dependence of the soliton depth, d, scaled to its ini- 
tial value, for T ~ 0.3 T c , v = 0.1c, and C12 = C22 = 0. This 
yields T\i% data points for the main graph (dashed lines) and 
the asymmetrical error bars obtained from 0.45 < d < 0.55 
(dotted lines). 



The main graph in Fig.|4]plots the temperature depen- 
dence of the half-lives for different initial soliton speeds. 
These reveal that the damping time decreases both with 
increasing temperature, as expected, and increasing ini- 
tial speed. The half-life is additionally found to decrease 
when collisions are taken into account, consistent with 
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the energy decay observed in Fig. [3] The latter change is 
dramatic for low temperatures and small initial velocities, 
although the effect of collisions becomes less important at 
higher temperatures. In particular, at high T and large 
initial v, the half- life is less than the period of the soliton 
oscillation (marked with a dashed line) even when colli- 
sions are absent, indicating that their inclusion will do 
little to further enhance the decay. Note that the large 
soliton damping observed here at relatively low temper- 
atures contrasts with the case of collective excitations, 
where finite T effects onl y b ecome important at higher 
temperatures, T > 0.5 T c [2l|]. 

This strong damping effect is present even though we 
have so far considered a rather idealized scenario with a 
large trap aspect ratio. In experiments the damping may 
be even more pronounced. To illustrate this, we explicitly 
consider the parameters of the experiment of Burger et 
al. 0: N = 1.5 x 10 5 87 Rb atoms, w z = 2ir x 14 Hz, 
u ± = 2tt x 425 Hz, and T ~ 0.5 T c 0]. 

Experimentally, a "phase-imprinting" method was 
used to generate solitons Q, where light was shone on 
one half of the condensate to create a phase imbalance, 
A<j>. To simulate this, we multiply our condensate wave- 
function at t = by , where <f>(z) — for z < -4/2, 
and 4>(z) — A(f> for z > 4/2. Experimentally, due to 
diffraction of the light, there is a smooth transition of 
the phase between the two halves, which takes place over 
a length of 4- We represent this transition with a linear 
ramp, (j> = (z/l e + 1/2) A<p, for -4/2 < z < 4/2. 

A further consideration when simulating the experi- 
ment is that the resulting solitons were too small to de- 
tect in situ, so the condensate was released from the trap 
and allowed to expand before imaging. We simulate this 
procedure by allowing the system to evolve in the trap for 
a variable time after phase imprinting, and subsequently 
setting the trap potential V = 0. The column density of 
the condensate is calculated t = 4 ms later, which simu- 
lates a typical absorption image obtained experimentally. 

Fig. [1] presents our computed expansion images for 
4 = 1 M m an d Acf) = it, for the cases of T = (top) and 
T ~ 0.5 T c (bottom). The corresponding experimental 
images for the first two columns (i.e. for times t < 12 ms) 
were presented in Ref. [7[ , where it was argued that the 
soliton had damped sufficiently to prevent it from being 
observed at subsequent times. On this timescale, our sim- 
ulated expansion images at finite T reveal that the soli- 
ton decays to approximately half of its original depth, in 
agreement with experiment. Moreover, in its subsequent 
evolution, the soliton spends a large amount of time near 
the edge of the trap where it is heavily damped. In fact, 
our finite temperature simulations reveal that the soliton 
never actually re-emerges from the edge of the trap, in 
stark contrast to the corresponding prediction of the GP 
equations (top row of Fig. [I}. Given the limited experi- 
mental visibility at the edge of the expanded profiles, our 
findings are thus in good agreement with the experiment 



of Ref. Q , as well as with the estimates of Ref. 18 1 



To conclude, we have studied the dynamics of dark soli- 
tons in elongated Bose Einstein condensates using finite 
temperature simulations. Unlike the undamped soliton 
oscillations predicted by the Gross-Pitaevskii equation, 
dissipation arises at finite temperatures, leading to an 
increase in the oscillation amplitude and a decrease in 
both the soliton depth and energy. This dissipation in- 
creases with temperature and initial velocity, and is in 
general sensitive to the inclusion of binary collisions be- 
tween the atoms. Nevertheless, our work indicates that 
the predicted dark soliton oscillations should be observ- 
able in realistic elongated geometries, provided that the 
temperature is a small fraction of the critical tempera- 
ture. Direct comparison of our simulations with the ex- 
periment of Burger et al. [7] shows very good agreement, 
demonstrating that thermal dissipation fully accounts for 
the absence of soliton oscillations in this experiment. 
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